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Perturbations to the homeostatic distribution of mechanical forces exerted by blood 
on the endothelial layer have been correlated with vascular pathologies including 
intracranial aneurysms and atherosclerosis. Recent computational work suggests 
that in order to correctly characterise such forces, the shear-thinning properties of 
blood must be taken into account. To the best of our knowledge, these findings have 
never been compared against experimentally observed pathological thresholds. In 
the current work, we apply the three-band diagram (TBD) analysis due to Gizzi et 
al. to assess the impact of the choice of blood rheology model on a computational 
model of the right middle cerebral artery. Our results show that the differences 
between the wall shear stress predicted by a Newtonian model and the well known 
Carreau-Yasuda generalized Newtonian model are only significant if the vascular 
pathology under study is associated with a pathological threshold in the range 
0.94 Pa to 1.56 Pa, where the results of the TBD analysis of the rheology models 
considered differs. Otherwise, we observe no significant differences. 

Keywords: Blood flow modelling, rheology, multiscale modelling, 
lattice-Boltzmann, three-band diagram analysis 

1. Introduction 

Physiology and medicine are being revolutionized by the growing role of information 
technology. Our ability to acquire and manage data on both animals and humans 
allows us to develop increasingly detailed computational models of the biological 
processes sustaining life. These models, together with the relevant experimental 
data, are helping researchers to gain insight into the physiology and pathology 
of the systems under study, in many cases beyond what is possible with purely 
observational methods. 

Cerebrovascular disorders, including intracranial aneurysms (ICA), which may 
rupture leading to subarachnoid haemorrhages, are one of the most prevalent and 
devastating diseases of adults, of worldwide concern. In the UK, the total bur- 
den has been estimated as £0.5 billion annually |Rivero- Arias et al., 2010| . It is 
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currently generally accepted that haemodynamics plays an important role in the 
appearance, evolution, and potential rupture of these type of vascular patholo- 
gies |Stehbens, 1989[ |Shojima et al., 2004) . More precisely, perturbations in the 
homeostatic distribution of mechanical forces exerted by the blood on the en- 
dothelial layer have been correlated not only to aneurysm initiation and rupture 
but also to the development of other vascular pathologies such as atherosclerosis 



Chatzizisis et al., 2008 



From a rheological point of view, blood is a shear-thinning fluid. This behaviour 
arises from the presence of red blood cells suspended in a medium known as blood 
plasma. Despite this fact, a large body of literature concerning computational 
haemodynamics characterises blood as a Newtonian fluid under the assumption 
that in large arteries the shear rate is large enough for the viscosity to be treated 



as effectively constant [Mejia et al., 2011 



There has been increasing interest during recent years in the comparison of 
blood rheology models using computational fluid dynamics (CFD) simulations in 
realistic computational domains reconstructed from medical image. These studies 
have been performed in both healthy vasculature (see e.g. Morbid ucci et al., 2011 



Box et al., 2005 Johnston et al., 2004|) and in the presence of aneurysms (see e.g 



Bernsdorf and Wang, 2009 Cavazzuti et al., 2011| |Fisher and Rossmann, 2009] ). 



In particular, authors pay special attention to the influence of the choice of rheology 
model on the computational estimates of wall shear stress (WSS) as a proxy for 
aneurysm rupture risk. It has been suggested that, in ICAs, the Newtonian simplifi- 
cation overestimates WSS (see e.g. Bernsdorf and Wang, 20 09. Xia ng et al., 2012] ) 



and may underestimate rupture risk. Furthermore, it has been argued that applica- 
tions requiring accurate WSS estimates (e.g. those concerning vascular remodelling 
and biomechanics) will suffer from modelling inaccuracy unless the generalized New- 
tonian (GN) properties of blood are taken into account. 

Different indices have been proposed as a way of integrating into CFD-based 
biomarkers of rupture risk the rich spatio-temporal structure of WSS induced by 
pulsatile flow in complex vascular networks. Minimal or maximal peak WSS, time- 
averaged WSS, and oscillatory shear index (OSI) |Ku et al., 1985| have been ex- 
tensively used, among others. In many of the studies cited previously, comparisons 
are performed based on colour map plots, of one or more of these indices, on the 
surface of the original computational domain. Differences between WSS estimates 
produced by different rheology models are hence presented in a very qualitative 
way. Furthermore, the link to actual rupture risk is based on the currently accepted 
low shear stress biomarker, but with no actual correlation to experimentally de- 
termined WSS thresholds. Several authors |Kallmes, 2012||Cebral and Meng, 2012] 
have recently criticised this approach to rupture risk quantification arguing that 
more quantitative methods are required in order to gain further insight into the 
problem. 

In a recent publication, Gizzi et al. |Gizzi et al., 20lT] proposed a new frame- 
work for the quantitative analysis of WSS: the so-called three-band diagram (TBD) 
analysis. The TBD analysis facilitates the evaluation of the WSS obtained at a 
given location over time (i.e. a WSS signal) against a range of WSS pathologi- 
cal thresholds (e.g. WSS magnitude lower than 0.5 Pa as reported in the case of 
atherosclerosis formation [Chatzizisis et al., 2008 ). The analysis determines how 
likely a given WSS signal is to be considered risky for any given threshold. Further- 
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Figure 1. 3D model of a subset of the right middle cerebral artery used in this work. Geom- 
etry segments are numbered A-E for later reference. The arrows indicate flow direction. 



more, the results of the analysis can be easily compared against pathological values 
of WSS observed experimentally. 

In the current work, we apply TBD analysis to assess the impact of the choice of 



blood rheology model on the WSS estimates of the HemeLB |Mazzeo and Coveney, 2008] 




lattice-Boltzmann blood flow solver in a high resolution three-dimensional model of 
the right middle cerebral artery. The rest of the paper is structured as follows: §2 
introduces the computational and mathematical models used in this work as well 
as the simulation workflow implemented, §3 presents the results of our simulation 
and their main implications, finally §4 summarises the main conclusions of the work 
and outlines future research directions. 



(i) Geometry generation 

The three-dimensional (3D) model of the middle cerebral artery used in this 
work (see figure [T]) is a subset of a geometrical model of the intracranial vascula- 
ture reconstructed from rotational angiography scans. It corresponds to a section 
of the right middle cerebral artery (MCA) in the vicinity of the internal carotid 
artery. The main geometrical features of the model are: i) vessels of variable di- 
ameter, ii) two bifurcations, and iii) vessel bending. The National Hospital for 
Neurology and Neurosurgery, London (UK) provided the original images in the 
framework of the GENIUS project [Mazzeo et al., 2010] as part of a larger dataset 
library. The dataset used in this work was segmented and the surface mesh in fig- 
ure [I] generated with the open source package Vascular Modeling Toolkit (VMTK) 
|Antiga et al., 2008| . 




2. Methods 



(a) 3D model of the middle cerebral artery 
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(ii) HemeLB 

HemeLB |Mazzeo and Coveney, 2008] is an open source software plat forrrjf] for 
modelling and simulation of blood flow in sparse vascular networks. It is comprised 
of tools for geometrical model preprocessing (i.e. regular grid volume meshing of 
surface meshes), simulation on massively parallel architectures, real-time visual- 
ization and steering, and data post-processing. To date, HemeLB has been suc- 
cessfully applied to the simulation of blood flow in healthy brain vasculature as 
well as in the presence of intracranial aneurysms. Particular attention has been 
paid to obtaining and presenting simulation results in a clinically meaningful way 
|Mazzeo et al., 20~10 . HemeLB uses the lattice-Boltzmann method for fluid dy- 



namics (see e.g. Chen and Doolen, 1998| ) since it allows efficient implementations 



in large-scale high performance computing infrastructures. For this work, we have 
developed an extension of HemeLB 's lattice Bhatnagar-Gross-Krook (LBGK) col- 
lision operator in order to accommodate both Newtonian and GN rheology models. 
We use the D3Q15 velocity set and the halfway bounce-back rule |Chen and Doolen, 1998] 



to enforce the no-slip boundary condition at the walls. We have recently shown Carver et al., 2012] 



that this combination of collision operator, velocity set, and wall boundary condi- 
tion performs well from both a numerical and a computational point of view in 
complex domains for typical blood flow Reynolds and Womersley numbers. 

(iii) Generalized Newtonian rheology 

The Carreau-Yasuda (CY) model is widely used to describe the shear-thinning 
behaviour of blood |Boyd et al., 2007, Ashrafizaadch and Bakhshaei, 2009 . In this 



model, the dynamic viscosity rj is related to the shear rate 7 through the following 
expression: 

Vti) =Voo + (m - ryoc) (1 + (A 7 ) a )^ , (2.1) 

where a, n, and A are empirically determined to fit a curve between regions of 
constant 77^ and rj . This model defines three different regimes: a Newtonian region 
with rjQ for low shear rate, followed by a shear-thinning region where r\ decreases 
with 7; finally, once rjao is reached a third Newtonian region with constant viscosity 
rjoo is defined for high shear rates. In this work, we will use the values provided 
in |Boyd et al, 2007| : 770 = 0.16 Pas, = 0.0035 Pas, A = 8.2 s, a = 0.64, and 
n = 0.2128 (both are dimensionless) . 

Figure [2] presents viscosity as a function of shear rate for the previous model and 
the Newtonian model considered in this work (77 = 3.5 x 10 _3 Pas). The Carreau- 
Yasuda model displays a smooth transition between ijo and 7700 . Significant haemo- 
dynamic differences between the two rheology models are expected for simulations 
with 7 < 100 s -1 . 

(b) ID model of the human vascular system 

In order to obtain inlet and outlet boundary conditions for our 3D model, we 
use a modified version of a one-dimensional (ID) model of the human vascular 
system previously published in Mu lder et al., 2011] . The original model includes 
a ID representation of the main arteries in the upper body (including the circle 



f The codebase is available under LGPL license from http://ccs.chem.ucl.ac.uk/hemelb 
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Figure 2. Dynamic viscosity 77 as a function of shear rate 7 for the Carreau-Yasuda (red) 

and the Newtonian (blue) models. 



Table 1. Geometrical characterisation of the 3D model used in this work. 



Segment 


Length 


Proximal end radius 


Distal end radius 


A 


8.9 mm 


1.30 mm 


1.25 mm 


B 


11.2mm 


1.25 mm 


1.10mm 


C 


13.6 mm 


1.10 mm 


0.85 mm 


D 


20.6 mm 


0.62 mm 


0.62 mm 


E 


21.1mm 


0.77 mm 


0.77 mm 



of Willis and both MCAs) and a zero-dimensional representation of the peripheral 
vasculature. The geometry is available as part of the open source software package 
pyNS Manini et al., 2012 which also implements a numerical solver for the ID 
pulse propagation mathematical model presented in |Huberts et al., 2012| . 

The original ID geometrical model does not include any detail about the vessels 
branching off the right MCA. Therefore, we modified it to include the two MCA 
side branches present in the 3D model (see figure [l}. In order to achieve this a 
ID characterisation of the 3D geometry was required. The open source software 
package VMTK was used to compute: a) the length of the centrelines associated 
with segments A-E in figure [II and b) the radii of the maximum inscribed spheres 
along the centrelines. Table [1] summarises the results and figure [3] shows a schematic 
of the ID model including the more detailed representation of the right MCA. 
The reader is referred to |Mulder et al., 20TT1 or the pyNS tool for the names and 
characteristics {e.g. length, radius, connectivity) of each of the segments of the 
vascular system. 
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Figure 3. ID model of the main arteries in the upper body. Each segment (in blue) rep- 
resents a different part of a human arterial system including characterisation of the 3D 
model of the right MCA in figure [T] (see top right arrow). Segment lengths are not to 
scale. The red dots represent: i) bifurcations, when found at the intersection of two or 
more segments, ii) zero-dimensional representations of the peripheral vasculature, when 
found at the end of an open-ended segment or iii) the heart in the case labelled with H. 



(c) Three-band diagram analysis 

The extraction of synthetic biomarkers of aneurysm rupture risk based on CFD 
analysis of patient-specific models is an active field of research. The three-band 
diagram analysis framework was recently proposed |Gizzi et al., 2011] as a way of 
generalising previously proposed indices (e.g. time-averaged WSS, OSI) and pro- 
viding a quantitative way of comparing temporal WSS signals against specific risk 
factors. 

Let t be the instantaneous traction vector at an arbitrary surface point with 
associated surface unit normal n such that 

ti{ Tjjiij , (2.2) 

where T is the deviatoric part of the full stress tensor of a fluid computed from 
the mesoscopic LB simulation variables as described in |Kruger et al., 2009| . The 
relationship between T and c, the full stress tensor, is 

= -PSij + Tij , (2.3) 

where P is the hydrodynamic pressure and 5 is the usual Kronecker delta tensor. 
Furthermore, for a generalized Newtonian incompressible fluid, T can be rewritten 
as 

T id = 27j(7)S y , (2.4) 
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where S is the shear rate tensor 



S ''i = \i d i v j + 9jVi) , (2.5) 



r\ is the dynamic viscosity of a fluid as a function of shear rate 7 = ^2SijSij, see 



e.g. equation (2.1 1, and v is the velocity vector. In the case of Newtonian fluids, 
vd) = V = const. 

In this work, we will consider the temporal evolution of the signed magnitude 
of the traction vector t (which we will refer as the WSS signal), i.e.: 

S(i):=s g n(t(t)-t)|t(t)|, (2.6) 

where t is the average traction vector over time, sgn is the sign function, and | • | is 
the magnitude of a given vector. 

Given a WSS signal, S(t), and a scalar risk factor a > 0, the TBD analysis 
defines a triplet of functions 

S+(t) := S(t)H(S{t) - a), 

S°(t) := S(t)H(a - S(t))H(a + S(t)), 

S-(t):=S(t)H(-S(t)-a), (2.7) 

where the Heaviside function H(x) is defined as H(x) = 1 if x > and H (x) = if 
x < 0. The closed support of the function S + is a set of time intervals of cardinality 
N + (a) (and similarly with S° and S~). The main idea behind the method is to 
inspect the number of such intervals as a function of the variable threshold a (i.e. 
the risk factor). It has been suggested that a WSS signal can be considered healthy if 
N + ' 0, ~(o-) > for a given risk factor a. The reader is referred to |Gizzi et al., 2011] 
for a more detailed description of the method. A Python implementation of the 
algorithm is freely available as part of HemeLB's postprocessing tools. 



(d) Simulation workflow 

Our simulation workflow is as follows. First, the original DICOM images were 
segmented, a region of interest was chosen, and a surface mesh was generated 
with the VMTK implementation of the marching cubes algorithm. A non-shrinking 
Taubin filter was then applied to smooth out imaging artifacts. Second, we loaded 
the resulting surface mesh in HemeLB's setup tool in order to choose the location 
of the inlet and the outlets and generate a regular grid discretisation of the volume 
enclosed by the surface (with a total of 4,161,046 grid points). Next, the pyNS 
solver was run to obtain pressure traces at inlet A and oulets C-E in figure [3j We 
then used these traces as inlet and outlet boundary conditions to run HemeLB 
simulations with both the Newtonian and Carreau-Yasuda rheology models for a 
total of three cardiac cycles. The following configuration paremeters were used in 
this work: timestep At — 2.5706 x 10 -6 s, spacestep Ax = 3.50 x 10~ 5 m, maxi- 
mum density difference in the domain Ap/p a — 0.019, Mach numbers (defined as 
the ratio between the largest velocity magnitude in the domain and lattice speed 
of sound c s = Jfa t ) 0.31 and 0.244 for the Newtonian and CY rheology models 
respectively, and LBGK relaxation parameter r = 0.522 for the Newtonian model 
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Table 2. Three-band diagram analysis sampling points. 
Point Coordinates 





(11.38, 36.89, 45.17) 


mm 


XI 


(10.98, 41.64, 47.97) 


mm 


X3 


(12.77, 41.15, 45.02) 


mm 


X4 


(12.34, 42.69, 46.17) 


mm 


X 5 


(11.88, 50.29, 41.31) 


mm 


X6 


(11.59, 50.55, 49.19) 


mm 


X T 


(12.89, 55.45, 44.75) 


mm 


Xs 


(12.82, 55.33, 45.58) 


mm 


Xg 


(14.57, 55.73, 52.63) 


mm 



and t = |+ Atc 2 p ^ or ^ c model (which becomes r = 0.522 for rjoo). Finally, we 
post-processed the simulation results to compute TBDs of the WSS signal at differ- 
ent geometrical locations over time. Table [2] shows their location for reproducibility 
purposes. 

The pyNS simulations were run on a single core of a 2 GHz Intel Core i7 laptop 
with 4GB of RAM and took on the order of minutes to run. The HemeLB simula- 
tions were run on 2048 cores of HECToR (UK's national supercomputer) and took 
53 min and 83 min for the Newtonian and Carreau-Yasuda models respectively. The 
difference in computing time between the two rheology models is due to the cost 



associated with computing 7 and evaluating equation (2.1 ) at each lattice site every 
timestep. The choice of core count is based on the number of lattice sites and the 
scalability analysis presented in |Groen et al., 2012b] in order to ensure efficient use 
of the computational resources. 

All the files required to run the pyNS and HemeLB simulations are available 
as part of the supplementary material associated with this paper. Both software 
packages are open source and freely available to the public. 

3. Results and discussion 

(a) Pressure profiles 

Figure 4] plots the pressure traces generated by pyNS at inlet A and outlets C-E 
of figure [3 It can be seen how the typical pressure variations observed throughout 
the cardiac cycle are well recovered. First, the sudden increase in pressure following 
the opening of the aortic valve after approximately 0.06 s, which peaks at around 
0.2 s and drops until approximately 0.4 s corresponds to the ventricular ejection 
phase. Secondly, the isovolumic relaxation phase happens between the previous 
time point and the mitral valve opening around 0.55 s. Finally, the ventricular 
filling phase occurs between the time the mitral valve opens and the end of the 
cycle at approximately 0.91 s (for a cardiac rate of roughly 66 beats per minute). 

(b) 3D simulations and three-band diagram analysis 

Figure [5] plots the traction vector t computed with HemeLB configured to use 
the CY rheology model at the end of diastole and at peak flow in systole (both in 
the third cardiac cycle simulated). The direction of t (which is consistent with the 
flow direction at any given point) and its magnitude follow a complex distribution 
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Figure 4. Pressure profiles at the inlet and outlets obtained from the ID model. Pressures 
are displayed relative to outlet D. Simulations were run for a total of three cardiac cycles. 

throughout the domain with notable changes during the cardiac cycle. The main 
characteristics are: i) areas with larger WSS can be found in the interior side of 
both branches leaving the bifurcation, ii) a zone of low WSS, or stagnation point, is 
found around the area where both branches meet, iii) there are substantial changes 
in the direction of t throughout the cardiac cycle around the stagnation point area, 
leading to oscillatory WSS signals. 

Figure [6] presents the TBD analysis at points 23, X4, and X5. The analysis con- 
firms the change of sign of the WSS signal (due to a change of direction in the 
traction field) occurring at X4 throughout the cardiac cycle (note the presence of 
a negative component in figures 6(c) and |6(d)j ). Furthermore, the negative compo- 
nent covers a wider range of threshold values in the Newtonian case. According to 
the criteria outlined in Giz zi et al., 2011] , an oscillatory WSS signal is considered 
healthy when the TBD analysis displays all three components above a given crit- 
ical threshold. In this case we observe a variation of just over 0.5 Pa between the 
predicted largest healthy threshold in the Newtonian and CY cases. 

In the case of X3 and x§ no negative component appears in the analysis indicating 
that there is no oscillation in the WSS signal. This rules out one of the main factors 
correlated with vascular disorders: oscillatory flow |Ku et al., 19 85*|. Nevertheless, 



such signals could still be considered risky if their mean value was below a given 



threshold value |Shojima et al., 2 004 1. We can clearly see in figures 6(a) 6(b) and 
6(e)| - 6(f) that the choice of rheology model would not have a significant impact on 



the assessment of risk since both diagrams are nearly identical. 

The points listed in table [2] which have not been analysed in figure [6] fall within 
one of the two previous scenarios. 



4. Conclusions 

In this work, we present a complete workflow for the simulation of blood flow in a 
patient-specific 3D model of the right middle cerebral artery. The main features of 
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(b) Peak flow during systole t = 1.97s. 



Figure 5. Traction vector t (estimated with the CY rheology model) at the upstream 
bifurcation of the 3D model. Vectors are scaled according to their magnitude. Points X3 
and X4 in table [2] are shown in green and magenta respectively. Visualizations generated 
with the open source software package Paraview Hender son, 2007] . 

the model are: i) the geometry is reconstructed from rotational angiography scans 
and discretised at high resolution (Ax = 3.5 x 10 -5 m), ii) inlet and outlet bound- 
ary conditions are obtained with a ID model of the complete vascular system, and 
iii) the rheological properties of the blood can be described with both Newtonian 
and generalized Newtonian models. Simulations run efficiently on the HECToR su- 
percomputer taking 53min to 83min for a domain comprising of 4,161,046 lattice 
sites. 

The workflow is applied to the comparison of two blood rheology models: a 
Newtonian model (77 = 3.5 x 10 _3 Pas) and the Carreau-Yasuda (CY) model. This 
is done in a quantitative manner in conjunction with the recently proposed three- 
band diagram (TBD) analysis framework. In agreement with previous work (see e.g. 
Morbiducci et al., 2011[ |Bernsdorf~a nd Wan g, 20*09] ), our results show variations 
in the haemodynamics recovered with each of the rheology models studied. However, 
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(a) Newtonian TBD at 23. 



(b) Carreau-Yasuda TBD at 23 





(c) Newtonian TBD at 24. 



(d) Carreau-Yasuda TBD at 24. 
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(e) Newtonian TBD at 25. 



(f) Carreau-Yasuda TBD at 25 



Figure 6. TBD analysis at points X3, £4, and x$. The values of '°' + ^(a) are 
presented as a stacked histogram. 



the evaluation of those results against a set of risk factors with TBD show little 
to no difference. In particular, a wall shear stress signal with strong oscillatory 
component was found to be risky for thresholds equal to or greater than 1.56 and 
0.94 Pa for the Newtonian and CY models respectively. In the case of non-oscillatory 
signals the analysis returns almost identical results in both cases. 

The main limitations of our study are: i) we have used a single geometry; in order 
to achieve statistical significance a larger number of other cases must be considered, 
and ii) the results presented may not hold in the presence of intracranial aneurysms 
(ICAs) or other vascular malformations where regions of much smaller shear rate 
may occur. 

Our future lines of work include applying the analysis methodology presented in 
the current work to a set of ICAs reconstructed from rotational angiography images. 
We also plan to investigate the coupling of the current 1D-3D haemodynamics 
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model with agent-based models of tissue remodelling using the multiscale modelling 
framework presented in |Groen et al., 2012a] . 
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